######################################################################## # Bivariate Distributions (APPL) Jeff Yang, Larry Leemis # dated 14 December 2008 #1111111111222222222233333333334444444444555555555566666666667777777777# # # Procedure name: None # # Other Procedures Called: None # # Date Revised: December 14, 2008 # # Purpose: Start up file # # Arguments: None # print(`PROCEDURES:`): print( `BiExpectedValue(XY, g(x, y))`, `BiTransform(XY, g, [h])`, `Covariance(XY)`, `Correlation(XY)`, `JointPDF(X, Y, [x], [y])`, `MarginalPDF(XY, x)`, `VerifyJointPDF(XY)`): print(``): print(`Procedure Notation:`): print(`XY is a bivariate random variable`): print(`Greek letters are numeric or symbolic parameters`): print(`x is numeric or symbolic`): print(`g and h are function`): print(`g(x, y) is an expression `): print(`Brackets [] denote optional parameters`): #1111111111222222222233333333334444444444555555555566666666667777777777# # # Procedure name: BiTransform(XY, g, [h]) # # Other Procedures Called: MarginalPDF # # Date Revised: December 14, 2008 # # Purpose: Calculates the joint PDF of UV via the bivariate # transformation-of-variables technique under the transformations # g and h; returns the marginal PDF of U # # Arguments: XY: continuous bivariate random variable # g: entered as a function (i.e. g := [(x,y)->x+y]) # h: same as g, default to h:= [(x,y)->y] # # with(VectorCalculus): BiTransform := proc(XY :: list(list), gp :: list, hp :: list) local g, h, sols, eqnx, eqny, constraints, PDF, i, UV,x,y,ncons,xinters,yinters,j,temp, geq, heq, u, v, totalconstraints, tempconsL, tempconsR, tempcons, endCons, endPDF, line1, line2, ib, jb, tempb, xcheck, ycheck, ucheck, vcheck, eqnxt, eqnyt,xcheck2, ycheck2,jac,p1,p2: g := gp; # # Check for the appropriate number of arguments # if (nargs <> 3 and nargs <> 2) then print(`ERROR(BiTransform): This procedure requires 2 or 3 arguments:`): print(`A bivariate random variable`): print(`A transformation`): print(`A dummy transformation (optional)`): RETURN(): fi: # # Check that the RV XY is in the list of 3 lists format # if (nops(XY) <> 3) then print(`ERROR(BiTransform): The RV XY must be a list of 3 sublists`): RETURN(): fi: # # Defaults dummy transformation to v = y # if(nargs=2) then h := [seq((x,y)->y,j=1..nops(XY[2]))]; # # Check that the number of transformations in g and h are the same # elif(nops(g) <> nops(hp)) then print(`ERROR(BiTransform): g and h must have the same number of partitions`): print(`equal to one or the number of partitions in XY`): RETURN(): else h:=hp; fi: # # apply g and h to all partitions # if((nops(g) = 1) and (nops(XY[2]) <> 1)) then g:=[seq(g[1],j=1..nops(XY[2]))]; h:=[seq(h[1],j=1..nops(XY[2]))]; fi: # # Check that the number of transformations matches the number of partitions # if((nops(g) <> 1) and (nops(g) <> nops(XY[2]))) then print(`ERROR(BiTransform): Number of constraints must equal number of transformations (unless given 1 transformation)`): RETURN(): fi: endPDF := []; endCons := NULL; for i from 1 to nops(XY[1]) do totalconstraints := []: x := 'x': y := 'y': xinters := []: # x intercept yinters := []: # corresponding y intercept line1 := []: # corresponding line 1 line2 := []: # corresponding line 2 ncons := nops(XY[2][i]): for j from 1 to ncons do temp := solve({lhs(XY[2][i][j])=rhs(XY[2][i][j]),lhs(XY[2][i][(j mod ncons)+1])=rhs(XY[2][i][(j mod ncons)+1])}); if(nops({temp}) > 1) then print(`ERROR(BiTransform)): Adjacent constraints intersect at two or more points`): RETURN(): elif(temp = NULL) then print(`ERROR(BiTransform): Adjacent constraints do not intersect`): RETURN(): fi: if((lhs(XY[2][i][j])=infinity or rhs(XY[2][i][j])=infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=infinity, y=0}: else temp := {x=0, y=infinity}: fi: elif((lhs(XY[2][i][j])=-infinity or rhs(XY[2][i][j])=-infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=-infinity, y=0}: else temp := {x=0, y=-infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=infinity or rhs(XY[2][i][(j mod ncons)+1])=infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=infinity, x=0}: else temp := {y=0, x=infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=-infinity or rhs(XY[2][i][(j mod ncons)+1])=-infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=-infinity, x=0}: else temp := {y=0, x=-infinity}: fi: fi: # # sets non "x,y" variables to x and y # if (lexorder(lhs(temp[1]),lhs(temp[2]))) then x:=lhs(temp[1]); y:=lhs(temp[2]); elif(lexorder(lhs(temp[2]),lhs(temp[1]))) then x:=lhs(temp[2]); y:=lhs(temp[1]); fi: if(temp <> NULL) then line1 := [op(line1), XY[2][i][j]]: line2 := [op(line2), XY[2][i][(j mod ncons)+1]]: fi: if(temp <> NULL and lhs(temp[1]) = x) then xinters := [op(xinters), rhs(temp[1])]: yinters := [op(yinters), rhs(temp[2])]: elif(temp <> NULL and lhs(temp[2]) = x) then xinters := [op(xinters), rhs(temp[2])]: yinters := [op(yinters), rhs(temp[1])]: fi: od: # #bubblesorts all 4 lists with respect to xinters # for ib to nops(xinters)-1 do for jb from ib+1 to nops(xinters) do if xinters[ib] > xinters[jb] then tempb := xinters[ib]: xinters[ib] := xinters[jb]; xinters[jb] := tempb: tempb := yinters[ib]: yinters[ib] := yinters[jb]; yinters[jb] := tempb: tempb := line1[ib]: line1[ib] := line1[jb]; line1[jb] := tempb: tempb := line2[ib]: line2[ib] := line2[jb]; line2[jb] := tempb: fi: od: od: # # find x and y values to check for correct inverse (assume .00001, ycheck is in support) # if(xinters[1] = -infinity) then xcheck := xinters[nops(xinters)]-0.0001; p1 := nops(xinters)-1; p2 := nops(xinters); else xcheck := xinters[1]+0.0001; p1 := 1; p2 := 2; fi: if(xinters[p1]=xinters[p2]) then ycheck := (yinters[p1]+yinters[p2])/2; else ycheck := (solve(subs(x=xcheck,lhs(line1[p1])=rhs(line2[p2])))+solve(subs(x=xcheck,lhs(line1[p1])=rhs(line2[p2])))) / 2 ; fi: ucheck := g[i](xcheck,ycheck); vcheck := h[i](xcheck,ycheck); geq := u = g[i](x,y); heq := v = h[i](x,y); sols:=allvalues(solve({geq,heq},{x,y})); if(sols = NULL) then print(`ERROR(BiTransform): The transformation equations have no solutions`): RETURN(): fi: # #Find correct inverse # if(nops([sols])=1) then sols := [sols]; fi: for j from 1 to nops([sols]) do #distinguish between x and y from solve if lhs(sols[j][1])=x then eqnxt := rhs(sols[j][1]): eqnyt := rhs(sols[j][2]): elif lhs(sols[j][2])=x then eqnxt := rhs(sols[j][2]): eqnyt := rhs(sols[j][1]): fi: #substitute back into x and y equations xcheck2:=unapply(eqnxt,[u,v])(ucheck,vcheck); ycheck2:=unapply(eqnyt,[u,v])(ucheck,vcheck); #if equal then update eqnx and eqny if(xcheck = xcheck2 and ycheck = ycheck2) then eqnx := eqnxt; eqny := eqnyt; fi: od: # #Remove redundant constraints # constraints := XY[2][i]; for j from 1 to nops(xinters) do if(g[i](xinters[j],yinters[j])=g[i](xinters[(j mod nops(xinters))+1],yinters[(j mod nops(xinters))+1]) and h[i](xinters[j],yinters[j])=h[i](xinters[(j mod nops(xinters))+1],yinters[(j mod nops(xinters))+1])) then if(j=nops(xinters)-1) then constraints := [op(constraints[1..nops(xinters)-1])]; elif(j=nops(xinters)) then constraints := [op(constraints[2..nops(xinters)])]; else constraints := [op(constraints[1..j]),op(constraints[j+2..nops(xinters)])]; fi: fi: od: # #cross multiply (assumes constraints are all inequalities) # for j from 1 to nops(constraints) do tempconsL := numer(lhs(subs({x=eqnx,y=eqny},constraints[j])))*denom(rhs(subs({x=eqnx,y=eqny},constraints[j]))); tempconsR := numer(rhs(subs({x=eqnx,y=eqny},constraints[j])))*denom(lhs(subs({x=eqnx,y=eqny},constraints[j]))); if(type(constraints[j],`<`) or type(constraints[j],`<=`) and ((tempconsL < 0 and tempconsL <0) or (tempconsL>0 and tempconsL>0))) then tempcons := tempconsL < tempconsR; else tempcons := tempconsL > tempconsR; fi: totalconstraints:= [op(totalconstraints),tempcons]; od: u:='u': v:='v'; assume(op(totalconstraints)); jac := Jacobian([eqnx,eqny],[u,v],'determinant')[2]; if(is(jac > 0)) then u:='u'; v:='v'; PDF := XY[1][i](eqnx,eqny) * Jacobian([eqnx,eqny],[u,v],'determinant')[2]: else u:='u'; v:='v'; PDF := XY[1][i](eqnx,eqny) * abs(Jacobian([eqnx,eqny],[u,v],'determinant')[2]): fi: endPDF := [op(endPDF),unapply(PDF,(u,v))]; endCons := endCons, totalconstraints; od: UV := [endPDF,[endCons],["Continuous","PDF"]]; #UV := [[unapply(PDF,(u,v))],[[op(totalconstraints)]],["Continuous","PDF"]]: return(MarginalPDF(UV,u)): end: #1111111111222222222233333333334444444444555555555566666666667777777777# # # Procedure name: BiExpectedValue(XY, g(x, y)) # # Other Procedures Called: None # # Date Revised: December 14, 2008 # # Purpose: Finds E[g(x, y)] # # Arguments: XY: continuous bivariate random variable # g(x, y): entered as an expression (i.e. a * x ^ 2 + 3 * y) # BiExpectedValue := proc(XY :: list(list), g) local NumLists, x, y, i, ncons, xinters, yinters, line1, line2, j, temp, ib, jb, tempb, absPDF, totalPDF, yupper, ylower, start, area, ind, lowerB, upperB; # # Check for the appropriate number of arguments # if (nargs <> 2) then print(`ERROR(BiExpectedValue): This procedure requires 2 arguments:`): print(`A bivariate random variable`): print(`A function g(x,y)`): RETURN(): fi: # # Check that the RV XY is in the list of 3 lists format # NumLists := nops(XY): if (NumLists <> 3) then print(`ERROR(BiExpectedValue): the RV XY must be a list of 3 sublists`): RETURN(): fi: # # Check that the third argument of XY is ["Continuous","PDF"] # if(XY[3] <> ["Continuous","PDF"]) then print(`ERROR(BiExpectedValue): Third argument of XY must be ["Continuous","PDF"]`): RETURN(): fi: totalPDF := 0: absPDF := 0: # #i loops through the number of segments of XY # for i from 1 to nops(XY[1]) do x := 'x': y := 'y': ncons := nops(XY[2][i]): xinters := []: # x intercept yinters := []: # corresponding y intercept line1 := []: # corresponding line 1 line2 := []: # corresponding line 2 # # j loops through the constraints for segment i # for j from 1 to nops(XY[2][i]) do temp := solve({lhs(XY[2][i][j])=rhs(XY[2][i][j]),lhs(XY[2][i][(j mod ncons)+1])=rhs(XY[2][i][(j mod ncons)+1])}); # # solve does not solve {x=0,y=+-infinity},{y=0,x=+-infinity} # if((lhs(XY[2][i][j])=infinity or rhs(XY[2][i][j])=infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=infinity, y=0}: else temp := {x=0, y=infinity}: fi: elif((lhs(XY[2][i][j])=-infinity or rhs(XY[2][i][j])=-infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=-infinity, y=0}: else temp := {x=0, y=-infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=infinity or rhs(XY[2][i][(j mod ncons)+1])=infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=infinity, x=0}: else temp := {y=0, x=infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=-infinity or rhs(XY[2][i][(j mod ncons)+1])=-infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=-infinity, x=0}: else temp := {y=0, x=-infinity}: fi: fi: if (nops({temp}) > 1) then print(`ERROR(BiExpectedValue): Adjacent constraints intersect at two or more points`): RETURN(): elif(temp = NULL) then print(`ERROR(BiExpectedValue): Adjacent constraints do not intersect`): RETURN(): fi: if (lexorder(lhs(temp[1]),lhs(temp[2]))) then x:=lhs(temp[1]); y:=lhs(temp[2]); elif(lexorder(lhs(temp[2]),lhs(temp[1]))) then x:=lhs(temp[2]); y:=lhs(temp[1]); fi: if (temp <> NULL) then line1 := [op(line1), XY[2][i][j]]: line2 := [op(line2), XY[2][i][(j mod ncons)+1]]: fi: if (temp <> NULL and lhs(temp[1]) = x) then xinters := [op(xinters), rhs(temp[1])]: yinters := [op(yinters), rhs(temp[2])]: elif(temp <> NULL and lhs(temp[1]) = y) then xinters := [op(xinters), rhs(temp[2])]: yinters := [op(yinters), rhs(temp[1])]: fi: if (nops(xinters) = ncons+1) then print(unbounded): fi: od: # #bubblesorts all 4 lists with respect to xinters # for ib to nops(xinters)-1 do for jb from ib+1 to nops(xinters) do if xinters[ib] > xinters[jb] then tempb := xinters[ib]: xinters[ib] := xinters[jb]; xinters[jb] := tempb: tempb := yinters[ib]: yinters[ib] := yinters[jb]; yinters[jb] := tempb: tempb := line1[ib]: line1[ib] := line1[jb]; line1[jb] := tempb: tempb := line2[ib]: line2[ib] := line2[jb]; line2[jb] := tempb: fi: od: od: # #default yupper and ylower assuming the start off is a vertical line # yupper := ({line1[1]} union {line2[1]}) intersect ({line1[2]} union {line2[2]}): ylower := ({line1[1]} union {line2[1]}) intersect ({line1[2]} union {line2[2]}): start := 1: # #start off from a point (figure out yupper and ylower) # if xinters[1] <> xinters[2] then area := int(solve(lhs(line1[1])=rhs(line1[1]),y)-solve(lhs(line2[1])=rhs(line2[1]),y), x=xinters[1]..xinters[2]): #PDF evaluated over the segment if(area > 0) then yupper := {line1[1]}: ylower := {line2[1]}: else ylower := {line1[1]}: yupper := {line2[1]}: fi: totalPDF := totalPDF + int(int(g*XY[1][i](x,y),y=solve(lhs(op(ylower))=rhs(op(ylower)),y).. solve(lhs(op(yupper))=rhs(op(yupper)),y),'continuous'), x=xinters[1]..xinters[2]): start := 2: #start position is two because first section has already been integrated fi: # #triangle case (left = point, right = line) # if(start = 2 and nops(xinters)=3 and xinters[3]=xinters[2]) then start:=4; fi: #print(`yupper,ylower,totalPDF,absPDF`,yupper,ylower,totalPDF,absPDF); # #begin calculating PDFs # for ind from start to nops(xinters) do #left xinters lie on a vertical line if xinters[ind] = xinters[ind+1] then #y ind < y ind+1 => ylower is other line intersecting ind and yupper is other line intersecting ind+1 if (yinters[ind] < yinters[ind+1]) then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: yupper := ({line1[ind+1]} union {line2[ind+1]}) minus yupper: else yupper := ({line1[ind]} union {line2[ind]}) minus yupper: ylower := ({line1[ind+1]} union {line2[ind+1]}) minus ylower: fi: # #infinity case # lowerB := solve(lhs(op(ylower))=rhs(op(ylower)),y): upperB := solve(lhs(op(yupper))=rhs(op(yupper)),y): if(lhs(op(yupper))=infinity or rhs(op(yupper))=infinity) then upperB := infinity: fi: if(lhs(op(ylower))=-infinity or rhs(op(ylower))=-infinity) then lowerB := -infinity: fi: totalPDF := totalPDF + int(int(g*XY[1][i](x,y),y=lowerB..upperB,'continuous'), x=xinters[ind]..xinters[ind+2]): ind := ind + 1: #left x is only 1 point elif xinters[ind] < xinters[ind+1] then #if ind+1 and ind don't lie on same upper/lower if op(({line1[ind]} union {line2[ind]}) intersect ({line1[ind+1]} union {line2[ind+1]})) = NULL then if yinters[ind] > yinters[ind+1] then yupper := ({line1[ind]} union {line2[ind]}) minus yupper: elif yinters[ind] < yinters[ind+1] then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: fi: #if ind+1 and ind do lie on same upper/lower else if op(yupper intersect ({line1[ind]} union {line2[ind]})) <> NULL then yupper := ({line1[ind]} union {line2[ind]}) minus yupper: elif op(ylower intersect ({line1[ind]} union {line2[ind]})) <> NULL then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: fi: fi: lowerB := solve(lhs(op(ylower))=rhs(op(ylower)),y): upperB := solve(lhs(op(yupper))=rhs(op(yupper)),y): if(lhs(op(yupper))=infinity or rhs(op(yupper))=infinity) then upperB := infinity: fi: if(lhs(op(ylower))=-infinity or rhs(op(ylower))=-infinity) then lowerB := -infinity: fi: totalPDF := totalPDF + int(int(g*XY[1][i](x,y),y=lowerB..upperB,'continuous'), x=xinters[ind]..xinters[ind+2]): fi: # #terminate if next is end or next two are equal and one of them is end # #print(`yupper,ylower,totalPDF,absPDF`,yupper,ylower,totalPDF,absPDF); if (ind+1 = nops(xinters)) or (xinters[ind+1]=xinters[ind+2] and ind+2=nops(xinters)) then break: fi: od: od: return(totalPDF): end: #1111111111222222222233333333334444444444555555555566666666667777777777# # # Procedure name: Correlation(XY) # # Other Procedures Called: MarginalPDF, Covariance, Variance (an APPL Procedure) # # Date Revised: December 14, 2008 # # Purpose: Finds the correlation of X and Y # # Arguments: XY: continuous bivariate random variable # Correlation := proc(XY :: list(list)) local temp,x,y; # # Check for the appropriate number of arguments # if (nargs <> 1) then print(`ERROR(Correlation): This procedure requires 1 argument:`): print(`A bivariate random variable`): RETURN(): fi: # # Check that the RV XY is in the list of 3 lists format # if (nops(XY) <> 3) then print(`ERROR(Correlation): the RV XY must be a list of 3 sublists`): RETURN(): fi: # # Check that the third argument of XY is ["Continuous","PDF"] # if(XY[3] <> ["Continuous","PDF"]) then print(`ERROR(Correlation): Third argument of XY must be ["Continuous","PDF"]`): RETURN(): fi: temp := solve({lhs(XY[2][1][1])=rhs(XY[2][1][1]),lhs(XY[2][1][2])=rhs(XY[2][1][2])}); x:=lhs(temp[1]); y:=lhs(temp[2]); return(Covariance(XY)/(Variance(MarginalPDF(XY,x))*Variance(MarginalPDF(XY,y)))^0.5); end: #1111111111222222222233333333334444444444555555555566666666667777777777# # # Procedure name: Covariance(XY) # # Other Procedures Called: BiExpectedValue, MarginalPDF, Mean (an APPL Procedure) # # Date Revised: December 14, 2008 # # Purpose: Finds the covariance of X and Y # # Arguments: XY: continuous bivariate random variable # Covariance := proc(XY :: list(list)) local temp,x,y; # # Check for the appropriate number of arguments # if (nargs <> 1) then print(`ERROR(Covariance): This procedure requires 1 argument:`): print(`A bivariate random variable`): RETURN(): fi: # # Check that the RV XY is in the list of 3 lists format # if (nops(XY) <> 3) then print(`ERROR(Covariance): The RV XY must be a list of 3 sublists`): RETURN(): fi: # # Check that the third argument of XY is ["Continuous","PDF"] # if(XY[3] <> ["Continuous","PDF"]) then print(`ERROR(Covariance): Third argument of XY must be ["Continuous","PDF"]`): RETURN(): fi: temp := solve({lhs(XY[2][1][1])=rhs(XY[2][1][1]),lhs(XY[2][1][2])=rhs(XY[2][1][2])}); x:=lhs(temp[1]); y:=lhs(temp[2]); RETURN(simplify(BiExpectedValue(XY,x*y)-Mean(MarginalPDF(XY,x))*Mean(MarginalPDF(XY,y)))): end: #1111111111222222222233333333334444444444555555555566666666667777777777# # # Procedure name: JointPDF(X, Y, [x], [y]) # # Other Procedures Called: None # # Date Revised: December 14, 2008 # # Purpose: Finds the marginal PDF of XY with respect to x # # Arguments: X and Y: continuous univariate random variable # x and y: symbolic variables for X and Y, respectively # JointPDF := proc(X :: list(list), Y :: list(list), xp, yp) local PDF, supp, i, j, x, y; # # Check for the appropriate number of arguments # if (nargs <> 4 and nargs <> 2) then print(`ERROR(JointPDF): This procedure requires 4 arguments:`): print(`Two univariate random variables`): print(`Two symbolic variables for the joint distribution (optional)`): RETURN(): fi: # # Defaults x and y to 'x' and 'y' # if(nargs = 4) then x := xp; y := yp; fi: # # Check that the RVs X and Y are in the list of 3 lists format # if (nops(X) <> 3 or nops(Y) <> 3) then print(`ERROR(JointPDF): The RVs X and Y must be in the list of 3 sublists format`): RETURN(): fi: # # Check that the RVs X and Y are PDFs # if (X[3] <> ["Continuous","PDF"] or Y[3] <> ["Continuous","PDF"]) then print(`ERROR(JointPDF): The third argument of X and Y must be ["Continuous","PDF"]`): RETURN(): fi: PDF := []; supp := []; for i from 1 to nops(X[1]) do for j from 1 to nops(Y[1]) do PDF := [op(PDF), unapply(simplify(X[1][i](x)*Y[1][j](y)),[x,y])]: supp := [op(supp), [x>X[2][i],y>Y[2][j],x 2) then print(`ERROR(MarginalPDF): This procedure requires 3 arguments:`): print(`A bivariate random variable`): print(`A marginal variable`): RETURN(): fi: # # Check that the RV XY is in the list of 3 lists format # NumLists := nops(XY): if (NumLists <> 3) then print(`ERROR(MarginalPDF): The RV XY must be a list of 3 sublists`): RETURN(): fi: # # Check that the third argument of XY is ["Continuous","PDF"] # if(XY[3] <> ["Continuous","PDF"]) then print(`ERROR(MarginalPDF): Third argument of XY must be ["Continuous","PDF"]`): RETURN(): fi: marginalPDF := []: suplower := []: supupper := []: # #i loops through the number of segments of XY # for i from 1 to nops(XY[1]) do x := 'x': y := 'y': ncons := nops(XY[2][i]): xinters := []: # x intercept yinters := []: # corresponding y intercept line1 := []: # corresponding line 1 line2 := []: # corresponding line 2 for j from 1 to nops(XY[2][i]) do temp := solve({lhs(XY[2][i][j])=rhs(XY[2][i][j]),lhs(XY[2][i][(j mod ncons)+1])=rhs(XY[2][i][(j mod ncons)+1])}); # # solve does not solve {x=0,y=+-infinity},{y=0,x=+-infinity} # if((lhs(XY[2][i][j])=infinity or rhs(XY[2][i][j])=infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=infinity, y=0}: else temp := {x=0, y=infinity}: fi: elif((lhs(XY[2][i][j])=-infinity or rhs(XY[2][i][j])=-infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=-infinity, y=0}: else temp := {x=0, y=-infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=infinity or rhs(XY[2][i][(j mod ncons)+1])=infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=infinity, x=0}: else temp := {y=0, x=infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=-infinity or rhs(XY[2][i][(j mod ncons)+1])=-infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=-infinity, x=0}: else temp := {y=0, x=-infinity}: fi: fi: if(nops({temp}) > 1) then print(`ERROR(MarginalPDF): Adjacent constraints intersect at two or more points`): RETURN(): elif(temp = NULL) then print(`ERROR(MarginalPDF): Adjacent constraints do not intersect`): RETURN(): fi: if(temp <> NULL) then line1 := [op(line1), XY[2][i][j]]: line2 := [op(line2), XY[2][i][(j mod ncons)+1]]: fi: if(temp <> NULL and lhs(temp[1]) = x) then y := lhs(temp[2]); xinters := [op(xinters), rhs(temp[1])]: yinters := [op(yinters), rhs(temp[2])]: elif(temp <> NULL and lhs(temp[2]) = x) then y := lhs(temp[1]); xinters := [op(xinters), rhs(temp[2])]: yinters := [op(yinters), rhs(temp[1])]: fi: od: # #bubblesorts all 4 lists with respect to xinters # for ib to nops(xinters)-1 do for jb from ib+1 to nops(xinters) do if xinters[ib] > xinters[jb] then tempb := xinters[ib]: xinters[ib] := xinters[jb]; xinters[jb] := tempb: tempb := yinters[ib]: yinters[ib] := yinters[jb]; yinters[jb] := tempb: tempb := line1[ib]: line1[ib] := line1[jb]; line1[jb] := tempb: tempb := line2[ib]: line2[ib] := line2[jb]; line2[jb] := tempb: fi: od: od: # #default yupper and ylower assume start off vertical line # yupper := ({line1[1]} union {line2[1]}) intersect ({line1[2]} union {line2[2]}): ylower := ({line1[1]} union {line2[1]}) intersect ({line1[2]} union {line2[2]}): start := 1: # #start off from a point (figure out yupper and ylower) # if xinters[1] <> xinters[2] then area := int(solve(lhs(line1[1])=rhs(line1[1]),y)-solve(lhs(line2[1])=rhs(line2[1]),y), x=xinters[1]..xinters[2]): #PDF evaluated over the segment if(area > 0) then yupper := {line1[1]}: ylower := {line2[1]}: else ylower := {line1[1]}: yupper := {line2[1]}: fi: marginalPDF := [op(marginalPDF),unapply(int(XY[1][i](x,y),y=solve(lhs(op(ylower))=rhs(op(ylower)),y).. solve(lhs(op(yupper))=rhs(op(yupper)),y),'continuous'),x)]: suplower :=[op(suplower),xinters[1]]: supupper :=[op(supupper),xinters[2]]: start := 2: fi: if(start = 2 and nops(xinters)=3 and xinters[2]=xinters[3]) then start:=4; fi: #print(`yupper,ylower,Marginal`,yupper,ylower,marginalPDF,`Range:`, xinters[1], xinters[2]); for ind from start to nops(xinters) do #left x start do not intersect (vertical line) if xinters[ind] = xinters[ind+1] then #y ind < y ind+1 => ylower is other line intersecting ind and yupper is other line intersecting ind+1 if(yinters[ind] < yinters[ind+1]) then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: yupper := ({line1[ind+1]} union {line2[ind+1]}) minus yupper: else yupper := ({line1[ind]} union {line2[ind]}) minus yupper: ylower := ({line1[ind+1]} union {line2[ind+1]}) minus ylower: fi: # #infinity case # marginalLB := solve(lhs(op(ylower))=rhs(op(ylower)),y): marginalUB := solve(lhs(op(yupper))=rhs(op(yupper)),y): if(lhs(op(yupper))=infinity or rhs(op(yupper))=infinity) then marginalUB := infinity: fi: if(lhs(op(ylower))=-infinity or rhs(op(ylower))=-infinity) then marginalLB := -infinity: fi: marginalPDF := [op(marginalPDF),unapply(int(XY[1][i](x,y),y=marginalLB..marginalUB ,'continuous'),x)]: ind := ind + 1: #left x only 1 point elif xinters[ind] < xinters[ind+1] then #if ind+1 and ind don't lie on same upper/lower if op(({line1[ind]} union {line2[ind]}) intersect ({line1[ind+1]} union {line2[ind+1]})) = NULL then if yinters[ind] > yinters[ind+1] then yupper := ({line1[ind]} union {line2[ind]}) minus yupper: elif yinters[ind] < yinters[ind+1] then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: fi: #if ind+1 and ind do lie on same upper/lower else if op(yupper intersect ({line1[ind]} union {line2[ind]})) <> NULL then yupper := ({line1[ind]} union {line2[ind]}) minus yupper: elif op(ylower intersect ({line1[ind]} union {line2[ind]})) <> NULL then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: fi: fi: marginalLB := solve(lhs(op(ylower))=rhs(op(ylower)),y): marginalUB := solve(lhs(op(yupper))=rhs(op(yupper)),y): if(lhs(op(yupper))=infinity or rhs(op(yupper))=infinity) then marginalUB := infinity: fi: if(lhs(op(ylower))=-infinity or rhs(op(ylower))=-infinity) then marginalLB := -infinity: fi: marginalPDF := [op(marginalPDF),unapply(int(XY[1][i](x,y),y=marginalLB..marginalUB ,'continuous'),x)]: fi: # #terminate if next is end or next two are equal and one of them is end # suplower :=[op(suplower),xinters[ind]]: supupper :=[op(supupper),xinters[ind+1]]: #print(`yupper,ylower,Marginal`,yupper,ylower,marginalPDF,`Range:`, xinters[ind], xinters[ind+1]); if (ind+1 = nops(xinters)) or (xinters[ind+1]=xinters[ind+2] and ind+2=nops(xinters)) then break: fi: od: od: # # Places individual parts in correct support # finalsupp := []: finalPDF := []: supp := sort([op(supupper),op(suplower)]): for i from 1 to nops(supp)-1 do if supp[i]<>supp[i+1] then finalsupp := [op(finalsupp),supp[i]]: finalPDF:= [op(finalPDF),unapply(0,x)]: for j from 1 to nops(marginalPDF) do if suplower[j] <= supp[i] and supupper[j] >= supp[i+1] then finalPDF[nops(finalPDF)]:= unapply(finalPDF[nops(finalPDF)](x)+marginalPDF[j](x),x): fi: od: fi: od: if finalsupp[nops(finalsupp)] <> supp[nops(supp)] then finalsupp:=[op(finalsupp),supp[nops(supp)]]: fi: # # Gets rid of type: piecewise # for i from 1 to nops(finalPDF) do finalPDF[i] := finalPDF[i](x) assuming (x>finalsupp[i] and x1 and finalPDF[nops(finalPDF)-1] = finalPDF[nops(finalPDF)]) then finalPDF := [op(finalPDF[1..nops(finalPDF)-1])]; finalsupp := [op(finalsupp[1..nops(finalPDF)-1]),op(finalsupp[nops(finalPDF)+1])]; fi: return([finalPDF,finalsupp, ["Continuous","PDF"]]); end: #1111111111222222222233333333334444444444555555555566666666667777777777# # # Procedure name: VerifyJointPDF(XY) # # Other Procedures Called: None # # Date Revised: December 14, 2008 # # Purpose: Verify that the PDF of XY is valid # # Arguments: XY: continuous bivariate random variable # VerifyJointPDF := proc(XY :: list(list)) local NumLists, x, y, i, ncons, xinters, yinters, line1, line2, j, temp, ib, jb, tempb, absPDF, totalPDF, yupper, ylower, start, area, ind, lowerB, upperB; # # Check for the appropriate number of arguments # if (nargs <> 1) then print(`ERROR(VerifyJointPDF): This procedure requires 1 argument:`): print(`A bivariate random variable`): RETURN(): fi: # # Check that the RV X is in the list of 3 lists format # NumLists := nops(XY): if (NumLists <> 3) then print(`ERROR(VerifyJointPDF): The RV XY must be a list of 3 sublists`): RETURN(): fi: # # Check that the third argument of XY is ["Continuous","PDF"] # if(XY[3] <> ["Continuous","PDF"]) then print(`ERROR(VerifyJointPDF): Third argument of XY must be ["Continuous","PDF"]`): RETURN(): fi: totalPDF := 0: absPDF := 0: # #i loops through the number of segments of XY # for i from 1 to nops(XY[1]) do x := 'x': y := 'y': ncons := nops(XY[2][i]): xinters := []: # x intercept yinters := []: # corresponding y intercept line1 := []: # corresponding line 1 line2 := []: # corresponding line 2 # # j loops through the constraints for segment i # for j from 1 to nops(XY[2][i]) do temp := solve({lhs(XY[2][i][j])=rhs(XY[2][i][j]),lhs(XY[2][i][(j mod ncons)+1])=rhs(XY[2][i][(j mod ncons)+1])}); # # solve does not solve {x=0,y=+-infinity},{y=0,x=+-infinity} # if((lhs(XY[2][i][j])=infinity or rhs(XY[2][i][j])=infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=infinity, y=0}: else temp := {x=0, y=infinity}: fi: elif((lhs(XY[2][i][j])=-infinity or rhs(XY[2][i][j])=-infinity) and (lhs(XY[2][i][(j mod ncons)+1])=0 or rhs(XY[2][i][(j mod ncons)+1])=0)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {x=-infinity, y=0}: else temp := {x=0, y=-infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=infinity or rhs(XY[2][i][(j mod ncons)+1])=infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=infinity, x=0}: else temp := {y=0, x=infinity}: fi: elif((lhs(XY[2][i][j])=0 or rhs(XY[2][i][j])=0) and (lhs(XY[2][i][(j mod ncons)+1])=-infinity or rhs(XY[2][i][(j mod ncons)+1])=-infinity)) then if(lhs(XY[2][i][j])=x or rhs(XY[2][i][j])=x) then temp := {y=-infinity, x=0}: else temp := {y=0, x=-infinity}: fi: fi: if (nops({temp}) > 1) then print(`ERROR(VerifyJointPDF): Adjacent constraints intersect at two or more points`): RETURN(): elif(temp = NULL) then print(`ERROR(VerifyJointPDF): Adjacent constraints do not intersect`): RETURN(): fi: if (lexorder(lhs(temp[1]),lhs(temp[2]))) then x:=lhs(temp[1]); y:=lhs(temp[2]); elif(lexorder(lhs(temp[2]),lhs(temp[1]))) then x:=lhs(temp[2]); y:=lhs(temp[1]); fi: if (temp <> NULL) then line1 := [op(line1), XY[2][i][j]]: line2 := [op(line2), XY[2][i][(j mod ncons)+1]]: fi: if (temp <> NULL and lhs(temp[1]) = x) then xinters := [op(xinters), rhs(temp[1])]: yinters := [op(yinters), rhs(temp[2])]: elif(temp <> NULL and lhs(temp[1]) = y) then xinters := [op(xinters), rhs(temp[2])]: yinters := [op(yinters), rhs(temp[1])]: fi: if (nops(xinters) = ncons+1) then print(unbounded): fi: od: # #bubblesorts all 4 lists with respect to xinters # for ib to nops(xinters)-1 do for jb from ib+1 to nops(xinters) do if xinters[ib] > xinters[jb] then tempb := xinters[ib]: xinters[ib] := xinters[jb]; xinters[jb] := tempb: tempb := yinters[ib]: yinters[ib] := yinters[jb]; yinters[jb] := tempb: tempb := line1[ib]: line1[ib] := line1[jb]; line1[jb] := tempb: tempb := line2[ib]: line2[ib] := line2[jb]; line2[jb] := tempb: fi: od: od: # #default yupper and ylower assuming the start off is a vertical line # yupper := ({line1[1]} union {line2[1]}) intersect ({line1[2]} union {line2[2]}): ylower := ({line1[1]} union {line2[1]}) intersect ({line1[2]} union {line2[2]}): start := 1: # #start off from a point (figure out yupper and ylower) # if xinters[1] <> xinters[2] then area := int(solve(lhs(line1[1])=rhs(line1[1]),y)-solve(lhs(line2[1])=rhs(line2[1]),y), x=xinters[1]..xinters[2]): #PDF evaluated over the segment if(area > 0) then yupper := {line1[1]}: ylower := {line2[1]}: else ylower := {line1[1]}: yupper := {line2[1]}: fi: totalPDF := totalPDF + int(int(XY[1][i](x,y),y=solve(lhs(op(ylower))=rhs(op(ylower)),y).. solve(lhs(op(yupper))=rhs(op(yupper)),y),'continuous'), x=xinters[1]..xinters[2]): absPDF := absPDF + int(int(abs(XY[1][i](x,y)),y=solve(lhs(op(ylower))=rhs(op(ylower)),y).. solve(lhs(op(yupper))=rhs(op(yupper)),y),'continuous'), x=xinters[1]..xinters[2]): start := 2: #start position is two because first section has already been integrated fi: # #triangle case (left = point, right = line) # if(start = 2 and nops(xinters)=3 and xinters[3]=xinters[2]) then start:=4; fi: #print(`yupper,ylower,totalPDF,absPDF`,yupper,ylower,totalPDF,absPDF); # #begin calculating PDFs # for ind from start to nops(xinters) do #left xinters lie on a vertical line if xinters[ind] = xinters[ind+1] then #y ind < y ind+1 => ylower is other line intersecting ind and yupper is other line intersecting ind+1 if (yinters[ind] < yinters[ind+1]) then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: yupper := ({line1[ind+1]} union {line2[ind+1]}) minus yupper: else yupper := ({line1[ind]} union {line2[ind]}) minus yupper: ylower := ({line1[ind+1]} union {line2[ind+1]}) minus ylower: fi: # #infinity case # lowerB := solve(lhs(op(ylower))=rhs(op(ylower)),y): upperB := solve(lhs(op(yupper))=rhs(op(yupper)),y): if(lhs(op(yupper))=infinity or rhs(op(yupper))=infinity) then upperB := infinity: fi: if(lhs(op(ylower))=-infinity or rhs(op(ylower))=-infinity) then lowerB := -infinity: fi: totalPDF := totalPDF + int(int(XY[1][i](x,y),y=lowerB..upperB,'continuous'), x=xinters[ind]..xinters[ind+2]): absPDF:= absPDF+int(int(abs(XY[1][i](x,y)),y=lowerB..upperB,'continuous'), x=xinters[ind]..xinters[ind+2]): ind := ind + 1: #left x is only 1 point elif xinters[ind] < xinters[ind+1] then #if ind+1 and ind don't lie on same upper/lower if op(({line1[ind]} union {line2[ind]}) intersect ({line1[ind+1]} union {line2[ind+1]})) = NULL then if yinters[ind] > yinters[ind+1] then yupper := ({line1[ind]} union {line2[ind]}) minus yupper: elif yinters[ind] < yinters[ind+1] then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: fi: #if ind+1 and ind do lie on same upper/lower else if op(yupper intersect ({line1[ind]} union {line2[ind]})) <> NULL then yupper := ({line1[ind]} union {line2[ind]}) minus yupper: elif op(ylower intersect ({line1[ind]} union {line2[ind]})) <> NULL then ylower := ({line1[ind]} union {line2[ind]}) minus ylower: fi: fi: lowerB := solve(lhs(op(ylower))=rhs(op(ylower)),y): upperB := solve(lhs(op(yupper))=rhs(op(yupper)),y): if(lhs(op(yupper))=infinity or rhs(op(yupper))=infinity) then upperB := infinity: fi: if(lhs(op(ylower))=-infinity or rhs(op(ylower))=-infinity) then lowerB := -infinity: fi: totalPDF := totalPDF + int(int(XY[1][i](x,y),y=lowerB..upperB,'continuous'), x=xinters[ind]..xinters[ind+2]): absPDF:= absPDF+int(int(abs(XY[1][i](x,y)),y=lowerB..upperB,'continuous'), x=xinters[ind]..xinters[ind+2]): fi: # #terminate if next is end or next two are equal and one of them is end # #print(`yupper,ylower,totalPDF,absPDF`,yupper,ylower,totalPDF,absPDF); if (ind+1 = nops(xinters)) or (xinters[ind+1]=xinters[ind+2] and ind+2=nops(xinters)) then break: fi: od: od: totalPDF := simplify(totalPDF): if(type(totalPDF, constant)=false) then print(`The area under f(x,y) evaluates to`, totalPDF): RETURN(): fi: if(evalf(totalPDF) < 0.9999999 or evalf(totalPDF) > 1.0000001) then print(`ERROR(VerifyJointPDF): The PDF of the given random variable`): print(`is NOT valid because the area under f(x,y) evaluates to`, totalPDF): RETURN(): elif ((evalf(totalPDF) > 0.9999999 and evalf(totalPDF) < 1.0000001) and (evalf(absPDF) > 0.9999999 and evalf(absPDF) < 1.0000001)) then print(`The area under f(x,y) evaluates to`, totalPDF): print(`f(x,y) is nonnegative`): RETURN(): else print(`ERROR(VerifyJointPDF): The PDF of the given random variable`): print(`is NOT valid because f(x,y) is negative for some value x,y`): print(`in its support`): RETURN(): fi: end: