weighted_statistics.dm42

export def wstat {
    if { MAT? } else {
        "X is not a matrix", AVIEW
        STOP
    }

    STO "A"
    GTO menu
}

def menu {
    CLMENU
    "mean", KEY 1 XEQ mean
    "std",  KEY 2 XEQ std
    "tsam", KEY 6 XEQ tsam
    FS? 03, KEY 9 GTO custom_menu
    
    MENU
    STOP
    GTO menu
}

def custom_menu {
    SF 27
    STOP
}

// Toggles flag 05
def tsam {
    if { FC? 05 } {
        SF 05
    } else {
        CF 05
    }
}

def index {
    INDEX "A"
    1, 1
    STOIJ
    DROPN 2
    J+, J- // For some reason, this is needed to clear flag 76
    // TODO: Assert N*2
}

// Sum(x * y) / Sum(y)
def mean {
    index()
    0, 0
    while { FC? 76} {
        RCLEL
        J+
        RCLEL
        STO+ ST T
        J-
        *
        +
        I+
    }
    X<>Y
    /
}

// Sqrt(Sum((x-mean)^2)*y))/Sum(y))
def std {
    mean()
    STO "B"
    DROP

    index()
    0, 0
    while { FC? 76} {
        RCLEL
        RCL "B"
        -
        X↑2
        J+
        RCLEL
        STO+ ST T
        J-
        *
        +
        I+
    }

    X<>Y
    if { FC? 05 } {
        1
        -
    }
    /
    SQRT
}