#OK to post homework
#Jingze Li, 2/3/2019,Assignment 3



#########3maple practice
> plot(sin(x),x=-2*Pi..2*Pi):
> plot([cos(t),1/2*sin(t),t=0..2*Pi]):
> plot([sqrt(x),3*log(x)],x=0..400):
> plot([sqrt(x),3*log(x)],x=0..400,numpoints=1000):
> L:=[[0,0],[1,1],[2,3],[3,2],[4,-2]]:
> plot(L):
> plot(L,style=point):
> p1:=plot([sqrt(x),3*log(x)],x=0..400,title=`the square root and log funcitons`):
> p2:=textplot([[360,16,`y=3log(x)`],[130,10,`y=sqrt(x)`]]):

> plot3d(exp(-(x^2+y^2-1)^2),x=-2..2,y=-2..2):
> x:='x'
                             x := x
;
> y
                               y
;
> y:='y'
                             y := y
;
> plot3d([sqrt(1+u^2)*cos(t),sqrt(1+u^2)*sin(t),u=-1..1,t=0..2*Pi]):
> sqrt(1+u^2)*cos(t),sqrt(1+u^2)*sin(t):
> with(plots):
> contourplot3d(exp(-(x^2+y^2-1)^5),x=-(1.3)..(1.3),y=-(1.3)..(1.3),filled=true,coloring=[blue,red]):
> with(plots):
> animate(1/(1+x*t),x=0..10,t=0..1,frames=100):
> animate3d([r*cos(t+a),r*sin(t+a),r^2*cos(2*t)+sin(a)],t=0..2*Pi,a=0..2*Pi,r=0..1,frames=10,style = patch,title=`the flying pizza`):
> v:=vector([1,2,3])
                         v := [1, 2, 3]
;
> print(v)
                           [1, 2, 3]
;
> A:=matrix(2,3,[a,b,c,d,e,f])
                              [a  b  c]                         A := [       ]                              [d  e  f]
;
> print(A)
                           [a  b  c]                           [       ]                           [d  e  f]
;
> op(A)
                           [a  b  c]                           [       ]                           [d  e  f]
;
> eval(A)
                           [a  b  c]                           [       ]                           [d  e  f]
;
> with(linalg):
> 
;


############solution to problem 3

#MomsAM(f,t,K): inputs a prob. gen. function and outputs the first K moments about the mean
 MomsAM:=proc(f,t,K) local L,i,F,ave:
 if subs(t=1,f)<>1 then
    print(`Not prob. dist.`):
    RETURN(FAIL):
 fi:
 
 F:=f:

  L:=[]:

   F:=t*diff(F,t):
    ave:=subs(t=1,F):

    F:=f/t^ave:

   for i from 1 to K do
    F:=t*diff(F,t):
    L:=[op(L),subs(t=1,F)]:
   od:
   [ave,op(2..nops(L),L)]:

 end:


#ScaledMomsAM(f,t,K): inputs a prob. gen. function and outputs the first K scaled moments about the mean
 ScaledMomsAM:=proc(f,t,K) local L,i,var,ave:
 if subs(t=1,f)<>1 then
    print(`Not prob. dist.`):
    RETURN(FAIL):
 fi:
 
  L:=[]:
  ave := MomsAM(f,t,1)[1]:
  var := MomsAM(f,t,2)[2]:
  
   for i from 3 to K do
    L:=[op(L),MomsAM(f,t,i)[i]/var^(K/2)]:
   od:
   [ave,var,op(L)]:

 end:





# this function computes the 1-14th moment of the prob-gen-func (1+t)/2)^n, as n takes different values.
#limS()  
limS:=proc(j) local i,n,L:

n:=[1, 10, 100, 1000, 10000,100000,200000,300000,400000,500000]:
#L:=matrix(14,nops(n)):
L:[]:

#for j from 1 to 14 do
for i from 1 to nops(n) do
  L:=[op(L),evalf(ScaledMomsAM(((1+t)/2)^n[i],t,j)[j])]:
od:
#od:

L:
end:
#####computation result:                                                                         guess limit point:
# j=1  .50000, 5., 50., 500., 5000., 50000., 100000., 150000., 200000., 250000                   inf
# j=2  .25000, 2.5000, 25., 250., 2500., 25000., 50000., 75000., 100000., 125000.                inf
# j=3  0., 0., 0., 0., 0., 0., 0., 0., 0., 0.                                                    0
# j=4  1., 2.8000, 2.9800, 2.9980, 2.9998, 3.0000, 3.0000, 3.0000, 3.0000, 3.0000                3
# j=5  0., 0., 0., 0., 0., 0., 0., 0., 0., 0.                                                    0
# j=6  1., 12.160, 14.702, 14.970, 14.997, 15.000, 15.000, 15.000, 15.000, 15.000                15
# j=7  0., 0., 0., 0., 0., 0., 0., 0., 0., 0.                                                    0
# j=8  1., 68.608, 100.86, 104.58, 104.96, 105.00, 105.00, 105.00, 105.00, 105.00                105
# j=9  0., 0., 0., 0., 0., 0., 0., 0., 0., 0.                                                    0
# j=10 1., 460.63, 883.62, 938.72, 944.37, 944.94, 944.97, 944.98, 944.98, 944.99                945
# j=11 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.                                                    0
# j=12 1., 3490.6, 9397.6, 10291., 10385., 10394., 10394., 10395., 10395., 10395.                10395
# j=13 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.                                                    0
# j=14 1., 28816., 117320., 133250., 134950., 135120., 135130., 135130., 135130., 135130.        135130 



##########solution to problem4

#this function competes the j'th moment of the random variable generated by quick sort
limQs:=proc(j) local i,L,n:


n:=[50,60,70,80,90,100]:

L:[]:


for i from 1 to nops(n) do
  L:=[op(L),evalf(ScaledMomsAM(PQs(n[i],t),t,j)[j])]:
od:

L:
end:

#when we let j=3, we get 1.0075, .99408, .98303, .97387, .96614, .95950
#when we let j=4, we get 4.4605, 4.4397, 4.4214, 4.4055, 4.3917, 4.3796
# these two sequences are both decreasing when n goes larger 
# thus, we need to comput more values to guess the limit point