examples/scripts/gradient.sh

Go to the documentation of this file.

#!/bin/zsh
script=${0##*/}

source $JPP_DIR/setenv.sh $JPP_DIR >& /dev/null
source $JPP_BIN/gradient.sh

set_variable:  DEBUG       SCRIPTS_DEBUG           3

if do_usage $*; then
    usage "$script [number of tests]"
fi

if (( $# == 1 )); then
    set_variable NUMBER_OF_TESTS      $1
else
    set_variable NUMBER_OF_TESTS      0
fi

#
# User function.
#
# \param 1     x
# \param 2     mean
# \param 3     sigma
# \param 4     signal
# \param 5     background
#
function gauss()
{
    let "PI = acos(-1.0)"
    let "u  = ($1 - $2) / $3"

    echo $(( $4 * exp(-0.5*u*u) / (sqrt(2.0*$PI) * $3)  +  $5 ))
}


#
# User data.
#
typeset -A    DATA                # (x,y) values


#
# Generate random data according true parameters.
#
function generate
{
    DATA=()

    for (( X = -5.0; $X <= +5.0; X += 0.5 )); do

        let "Y  = `gauss $X $TRUE_MEAN $TRUE_SIGMA $TRUE_SIGNAL $TRUE_BACKGROUND`"
        
        DATA[$X]=`$JPP_DIR/examples/JMath/getPoisson -e $Y`
    done
}


#
# Get chi2.
#
# The function is evaluated using parameters MEAN, SIGMA, SIGNAL and BACKGROUND.
#
# \param 1     variable containing chi2 value on return
#
function getChi2()
{
    let "SUM = 0.0"

    for X Y in ${(kv)DATA}; do

        let "F  = `gauss $X $MEAN $SIGMA $SIGNAL $BACKGROUND`"
        let "DY = ($F - $Y) * ($F - $Y)"

        if (( $Y > 0 )); then
            let "DY /= $Y"
        else
            let "DY  = $F"
        fi

        let "SUM += $DY"

    done

    eval $1=$SUM
}


# True values.

let "TRUE_MEAN        =    0.0"
let "TRUE_SIGMA       =    1.0"
let "TRUE_SIGNAL      = 1000.0"
let "TRUE_BACKGROUND  =   10.0"


#
# Print parameters.
#
function typeout()
{
    for KEY in MEAN SIGMA SIGNAL BACKGROUND; do

        set_variable P1 $KEY
        set_variable P2 TRUE_$KEY
        
        printf "%-12s  %12.3f (fit) %12.3f (true)\n"  $KEY  ${(P)P1}  ${(P)P2}
    done
}


if (( 1 )); then

    DATA=(
        -3.00000     16.00000
        -2.50000     22.00000
        -2.00000     70.00000
        -1.50000    152.00000
        -1.00000    261.00000
        -0.50000    337.00000
         0.00000    378.00000
         0.50000    360.00000
         1.00000    253.00000
         1.50000    129.00000
         2.00000     72.00000
         2.50000     22.00000
         3.00000     11.00000)

    # Fit parameters.

    PARAMETERS=()
    
    PARAMETERS['let "MEAN        +=  %"']=1.0e-2
    PARAMETERS['let "SIGMA       +=  %"']=1.0e-2
    PARAMETERS['let "SIGNAL      +=  %"']=5.0e-0
    PARAMETERS['let "BACKGROUND  +=  %"']=5.0e-1

    # Start values.

    let "MEAN        =    0.5"
    let "SIGMA       =    0.5"
    let "SIGNAL      =  700.0"
    let "BACKGROUND  =    0.0"

    gradient 0  Y1

    printf "chi2/NDF      %12.5f (fit) (gradient)\n" $(($Y1 / (${#DATA} - ${#PARAMETERS})))

    simplex  Y1

    printf "chi2/NDF      %12.5f (fit) (simplex) \n" $(($Y1 / (${#DATA} - ${#PARAMETERS})))

    typeout

    let "MEAN        = $TRUE_MEAN"
    let "SIGMA       = $TRUE_SIGMA"
    let "SIGNAL      = $TRUE_SIGNAL"
    let "BACKGROUND  = $TRUE_BACKGROUND"

    getChi2  Y2

    printf "chi2/NDF      %12.5f (fit) %12.5f (true)\n" $(($Y1 / (${#DATA} - ${#PARAMETERS})))  $(($Y2 / ${#DATA}))
fi


if (( $NUMBER_OF_TESTS > 0 )); then

    let "CHI2_FIT  = 0.0"
    let "CHI2_TRUE = 0.0"

    for (( m = 0; $m != $NUMBER_OF_TESTS; ++m )); do

        generate

        # Fit parameters.

        PARAMETERS=()
    
        PARAMETERS['let "MEAN        +=  %"']=1.0e-2
        PARAMETERS['let "SIGMA       +=  %"']=1.0e-2
        PARAMETERS['let "SIGNAL      +=  %"']=5.0e-0
        PARAMETERS['let "BACKGROUND  +=  %"']=5.0e-1

        # Start values.

        let "MEAN        =    0.5"
        let "SIGMA       =    0.5"
        let "SIGNAL      =  700.0"
        let "BACKGROUND  =    0.0"

        gradient 0  Y1

        let "MEAN        = $TRUE_MEAN"
        let "SIGMA       = $TRUE_SIGMA"
        let "SIGNAL      = $TRUE_SIGNAL"
        let "BACKGROUND  = $TRUE_BACKGROUND"
        
        getChi2  Y2
        
        printf "chi2/NDF  %3d %12.5f (fit) %12.5f (true)\n" $m  $(($Y1 / (${#DATA} - ${#PARAMETERS})))  $(($Y2 / ${#DATA}))

        let "CHI2_FIT  += $Y1 / (${#DATA} + ${#PARAMETERS})"
        let "CHI2_TRUE += $Y2 /  ${#DATA}"
    done

    printf "<>            %12.5f (fit) %12.5f (true)\n"  $(($CHI2_FIT / $NUMBER_OF_TESTS))  $(($CHI2_TRUE / $NUMBER_OF_TESTS))
fi