Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to express a simple probability problem in Mathematica, but am having trouble getting my calculations to execute at a reasonable speed.

I have an object whose unknown location is modeled as a random variable with coordinates $\Theta=(\Theta_1,\Theta_2)$ distributed as

$$f_\Theta(\theta)=f_\Theta(\theta_1,\theta_2)=\frac{1}{2\pi}e^{-(\theta_1^2+\theta_2^2)/2}$$

and a collection of sensors at given locations $(a_i, b_i)$ each of which has a probability of detecting a signal from the object that decays exponentially with the distance, modeled as a boolean-valued random variable $X_i$ with

$$\mathbf{P}(X_i=T|\Theta=(\theta_1,\theta_2))=p_{X_i|\Theta}(T|\theta)=e^{-\sqrt{(a_i-\theta_1)^2+(b_i-\theta_2)^2}}$$

so that the joint PDF of a set of sensor states $x=\{x_1,x_2,\ldots\}$ and object location $\theta$ is

$$f_{X,\Theta}(x,\theta)=f_\Theta(\theta)p_{X|\Theta}(x|\theta)=f_\Theta(\theta)\prod_{i|x_i=T}p_{X_i|\Theta}(T|\theta)\prod_{i|x_i=F}(1-p_{X_i|\Theta}(T|\theta))$$

so that

$$p_X(x)=\int_{-\infty}^\infty\int_{-\infty}^\infty f_{X,\Theta}(x,\theta)d\theta_1 d\theta_2$$

and thus the probability that the object is at location $(\theta_1,\theta_2)$ given a set of sensor observations $x$ is given by the PDF

$$f_{\Theta|X}(\theta|x)=\frac{f_{X,\Theta}(x,\theta)}{p_X(x)}$$

I've tried to express this in Mathematica as

PDetect[sensloc_, srcloc_] := Exp[-EuclideanDistance[sensloc, srcloc]];
PDFSrcloc[srcloc_] := 1/(2 Pi)Exp[-(srcloc[[1]]^2 + srcloc[[2]]^2)/2];
                      (*PDF[BinormalDistribution[{1,1},0],srcloc]*)
PDFJointSrcSens[srcloc_, senslocs_, sensstates_] := 
  PDFSrcloc[srcloc]*
   If[Length[sensstates] > 0, 
    Times @@ MapThread[
      If[#1, PDetect[#2, srcloc],  1 - PDetect[#2, srcloc]] &, {sensstates, senslocs}], 1];
PDFSensLocs[senslocs_, sensstates_] := 
  PDFSensLocs[senslocs, sensstates] = 
   If[Length[sensstates] > 0, 
    NIntegrate[
     PDFJointSrcSens[{x, y}, senslocs, sensstates], {x, -Infinity, 
      Infinity}, {y, -Infinity, Infinity}], 1];
PDFSrcLocGivenSens[srcloc_, senslocs_, sensstates_] := (PDFJointSrcSens[srcloc, senslocs, sensstates])/
   PDFSensLocs[senslocs, sensstates];

This seems to give me the results I expect, but performance is quite slow and plagued with beeps and warnings for many values. This is adequate for simple evaluations and marginal plots, but I need to put this in a Manipulate in which I dynamically plot the DensityPlot for PDFSrcLocGivenSens — that is, $f_{\Theta|X}(\theta|x)$ — as I add and remove sensors (adding/removing and changing the values of $x_i$; adding and removing elements to senslocs and adjusting the corresponding sensstates).

What can I do to speed up my calculations and avoid the beeps I'm getting?


I'm (clearly) not great at this, so any assistance in a more elegant formulation would be appreciated. I'm particularly unhappy with my argument and function names.

The original formulation of the scenario, which I may have garbled, can be found on p. 418 of Introduction to Probability, 2nd Edition, by Dimitri P. Bertsekas and John N. Tsitsiklis.

share|improve this question
1  
In addition to the source code, could you post (or provide a link to) sample values for sensloc and srcloc? That would allow us to see what the performance looks like and get an idea of what kind of error messages you're receiving. If it's too large to fit in your question post comfortably, I recommend posting it at pastebin.com . –  Pillsy Apr 12 at 14:38
    
@Pillsy: Here's the code in context, where you can generate typical values for each. In general, srcloc will simply be a point such as {1.4, -0.5}, senslocs will be a list of such points, and sensstates will be a boolean list of the same length. –  raxacoricofallapatorius Apr 12 at 15:04

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.