vka725 has asked for the wisdom of the Perl Monks concerning the following question:

Hi everyone!

I'm currently working on writing a little script to be able to calculate the average energy and the standard deviation of that energy of these molecular simulations I ran. There's about 600 snapshots for each molecule and I have to find the average and standard deviation in the energy of the snap shots for each molecule.

I am not only a super Perl newb, i'm new to coding as well. My PI helped me out for the first part (using it to find averages) and it worked. The second part was adding to it to be able to calculate the standard deviation.

So I was able to obtain some kind of standard deviation for 3 snapshots but they did not match the standard deviation i used in excel (i used STDEV.P since we have the entire population, N, not N-1).

I can't seem to get it right and I want to see if anyone can help me out here before I go and ask my PI. I want to understand the logic behind this as well, rather than have him just type something out for me.

Thank you! 
#!/usr/bin/env perl

$conf = $ARGV[0];
$n1 = $ARGV[1];
$n2 = $ARGV[2];

$esum0 = 0;
$esum1 = 0;
for($i=$n1;$i<=$n2;$i++){

    $i = sprintf("%0.3d",$i);


    system("analyze $conf.$i E | grep \"ACE\"  -A 1 > log ");


    open(IN,"log");
    @temp = split(" ",<IN>);
    $e0 = $temp[2];
    @temp = split(" ",<IN>);
    $e1 = $temp[1];
    close IN;



   print "$i  $e0 $e1 
";


#---average-calculation---

   $esum0 = $esum0 + $e0;
   $esum1 = $esum1 + $e1;

}

{
   $esum0=$esum0/$n2;
   $esum1=$esum1/$n2;
    $esum0= sprintf("%0.4f",$esum0);
   $esum1= sprintf("%0.4f",$esum1);

}

#averages
my@coolbeans = ('averages',  " $esum0", "$esum1");

print "
@coolbeans\n";




#---standard-deviation---


for($i=$n1;$i<=$n2;$i++){

        $estd0=($esum0-$e0) ** 2;
        $estd0=($estd0/$n2) **0.5;
        $estd0= sprintf ("%0.4f", $estd0);

        $estd1=($esum1-$e1) ** 2;
        $estd1=($estd1/$n2) **0.5;
        $estd1= sprintf ("%0.4f", $estd1);
}

#stdev
my@arroba = ('stdev', "$estd0", "$estd1");
print "
@arroba\n";

[download]

Replies are listed 'Best First'.
Re: standard deviation help
by GotToBTru (Prior) on Jul 07, 2015 at 19:56 UTC
[reply]
      thank you!
[reply]
Re: standard deviation help
by Laurent_R (Canon) on Jul 07, 2015 at 21:02 UTC
    For fun rather than for serious purpose, using the reduce function of the List::Utils module (or a home-made reduce function*) to build a library of statistical utilities in functional programming style. 
    sub sum {
    reduce { $a + $b } 0, @_; 
}
sub avg {     
    return undef unless @_; # avoid division by 0
    sum (@_)/scalar @_;
}
sub variance {
    return undef unless @_; # avoid division by 0
    my $avg = avg (@_);
    sum ( map {($_ - $avg)**2} @_ )/ scalar @_;
}
sub std_dev {   
    sqrt variance (@_);
}

    [download]
    The reduce utility can also be written in pure Perl: 
    sub reduce (&@) {
    my $code_ref = shift;
    my $result = shift;
    for (@_ ) {
        local ($a, $b) = ($result, $_ );
        $result = $code_ref->($a, $b )
    }
    return $result;
}

    [download]
[reply]
[d/l]
[select]
Re: standard deviation help
by hippo (Archbishop) on Jul 07, 2015 at 21:39 UTC
    I am not only a super Perl newb, i'm new to coding as well.

    Since there are two good replies to your specific questions already, permit me to digress and give you the best and most concise piece of general advice I can.

    Always use strict. Always.

    If you just add that one pragma right below the #! your code will no longer compile and you'll have to go through and make a bunch of changes so that it will compile. That may seem like pointless housekeeping but as someone who is not so new to perl or to programming let me impress upon you that this will save your bacon over and over and over again.

[reply]
      perfect! thank you so much for the tip :-)
[reply]
Re: standard deviation help
by poj (Abbot) on Jul 08, 2015 at 10:53 UTC
    ..the first part (using it to find averages) and it worked.

    Presumably that was with $n1 = 1 otherwise I can't see how $esum0=$esum0/$n2 would work.

    Shouldn't it be $eavg0 = $esum0/($n2-$n1+1); ?

    poj
[reply]
[d/l]
[select]

Back to Seekers of Perl Wisdom