number_simulations := 100000:
profit_array := Array(1..number_simulations, 1..6):
for kits_ordered from 100 to 225 by 25 do
j := (kits_ordered - 100)/25 + 1;
for i from 1 to number_simulations by 1 do
u := rand()/10^12;
if (u<.2) then kits_demand := 100 elif (u<.4) then kits_demand := 125 elif (u<.7) then kits_demand := 150 elif (u<.9) then kits_demand := 175 elif (u<.95) then kits_demand := 200 else kits_demand := 225 end if;
if (kits_ordered < kits_demand) then profit_array[i,j] := 10*kits_ordered else profit_array[i,j]:= (10*kits_demand - 10*(kits_ordered - kits_demand)) end if;
end do:
end do:
averages := Array(1..6):
for j from 1 to 6 by 1 do
averages[j] := evalf(add(profit_array[i,j], i=1..number_simulations)/number_simulations);
end do:
averages;
standard_deviations := Array(1..6):
for j from 1 to 6 by 1 do
standard_deviations[j] := sqrt(evalf(add( (profit_array[i,j] - averages[j])^2, i=1..number_simulations)/(number_simulations -1)));
end do:
standard_deviations;