; --------------- STATISTICS ----------------
; regression Y=A.F(x,T)+B with T optimization
; -------- guillaume.tello@orange.fr --------

; procedures

#label A  edit_values
#label B  set_points_number
#label C  get_correlation
#label E  function
#label A' init_T_param
#label B' set_precision
#label C' optimize_T

; variables

				; REG01 to REG06 for statistics
#reg  7 I                       ; pointer for "points" array
#reg  8 precision               ; when optimizing T
#reg  9 counter                 ; loop when entering stats in memory
#reg 10 T                       ; parameter to optimize
#reg 11 delta                   ; value 1.0001 to approximize f'() and f"()
#reg 12 f"                      ; f"()=Correl(deltaT)-Correl(delta.T)
#reg 13 f'                      ; f'()=Correl(delta.T)-Correl(T)
#reg 14 points_nbr              ; number of x;y points in "points" array

; arrays

#array points 15 99             ; array starting at 15 with no limit (except memory)

; programs

#pgm 1 master_diagnostic        ; this program...
#label CLR stat_init            ; ... to clear statistic registers (1 to 6)

;-----------------
; A to edit values
;-----------------

LBL edit_values

    BOT points STO I            ; pointer to the start of the Array
    DO
        RCL points(I) R/S       ; and set/edit every value
        STO points(I) 
        INC I
    LOOP                        ; no limit

;-------------------------------------
; C to get the correlation coefficient
;-------------------------------------

LBL get_correlation

    PGM master_diagnostic       ; init stats registers
        SBR stat_init
        
    I={ BOT points }            ; start of array
    counter={ @points_nbr }     ; to loop
    
    FOR counter
    
        @points(I) SBR f(x)     ; computes f(x) for each x
        X/T
        INC I 
        @points(I) S+           ; and add (f(Xi);Yi) as a new value
        INC I
        
    LOOP
    
    Correlation                 ; returns correlation coefficient (OP 13)
    
RTN

;------------------------------------------------------------
; optimize T to get the higher correlation coefficient
; this is done considering the coefficient as a function of T
; and looking for d.coeff/d.T = 0 with a Newton approximation
;------------------------------------------------------------

LBL optimize_T

    DO
    
        T*{ @delta X }                         ; T=1.0001.T
            get_correlation STO f"              ; f"()=Correl(delta.T)
        
        T/{ @delta }                            ; T=1.0001.T
            get_correlation INV SUM f"          ; f"()=Correl(delta.T)-Correl(delta.T)
                            STO f'              ; f'()=Correl(delta.T)
                             
        T/{ @delta }                            ; T=T   
            get_correlation INV SUM f'          ; f'()=Correl(delta.T)-Correl(T)
        
        @precision X/T                          ; for comparison
    
        T-{ @f" / @f' - 1 = 1/X                 ; computes dx=.0001T/(f"/f'-1) (newton for f'=0)
            * @T * $-4 10^x = }                 ; and then T=T-dx
        ABS                                     ; and ABS(dx) as precision met
        
    UNTIL(x<t)                                  ; loop until ABS(dx)<precision
    
    Correlation                                 ; Correlation in T
    X/T
    RCL T                                       ; and best parameter in X

RTN

;---------------------------------------------
; T A' to set the initial value of T parameter
;---------------------------------------------

LBL init_T_param

    STO T

RTN

;---------------------------------------
; e B' to set precision (and init delta)
;---------------------------------------

LBL set_precision

    X/T
    delta={ $1.0001 }
    X/T STO precision

RTN

;---------------------------
; x D to get Y for a given x
;---------------------------

LBL D

    function Y(x)

RTN

;--------------------------------
; n B to set the number of points
;--------------------------------

LBL set_points_number

    STO points_nbr

RTN

;-----------------------------------------
; user function, X at the display, t in @T
;        Reg00 free for the user if needed
;-----------------------------------------

LBL function

f(x):
    X/T @T Y^X X/T =

RTN

#end
