Input LogLikelihoodFunction directly

[spasetto]

Hi, Is it possible to input directly the LogLikelihoodFunction in defineInferenceProblem? It seems to be in the options (just looking at defineInferenceProblem[]) but I'm struggling with the syntax: From your example code:

BlockRandom[ SeedRandom[1]; testdata = (# -> (# + 2 RandomReal[])) & /@ RandomReal[{-1, 1}, 20]; ListPlot[List @@@ testdata] ]

instead of

obj = defineInferenceProblem[ "Data" -> testdata, "GeneratingDistribution" -> NormalDistribution[ a x + b, sigma ], "Parameters" -> {{a, -10, 10}, {b, -10, 10}, {sigma, 0.1, 10}}, "IndependentVariables" -> {x}, "PriorDistribution" -> {"LocationParameter", "LocationParameter", "ScaleParameter"} ]

that in result = nestedSampling[obj] gives me something like a=0.90... b=0.88..., I would need to be able to write something like:

obj = defineInferenceProblem[ "Data" -> testdata, "LogLikelihoodFunction" -> LogLikelihood[ NormalDistribution[0, sigma], Table[testdata[[i, 2]] - (a testdata[[i, 1]] + b), {i, 1, Length[testdata]}]], "Parameters" -> {{a, -10, 10}, {b, -10, 10}, {sigma, 0.1, 10}}, "PriorDistribution" -> {"LocationParameter", "LocationParameter", "ScaleParameter"} ]

...where at a certain point also the function ax+b should be entered(?). But I always failed, inferenceObject[$Failed].

Thx for your help!

PS: I've tried also something with directLogLikelihoodFunction, but the results have no sense.

[SjoerdSmitWolfram]

Hi,

I can see what's happening here. The "LogLikelihoodFunction" you're giving is not actually a function; it's an expression. A function is something that takes arguments and returns a value. It should work if you use:

"LogLikelihoodFunction" -> expressionToFunction[
 LogLikelihood[NormalDistribution[0, sigma], Table[testdata[[i, 2]] - (a testdata[[i, 1]] + b), {i, 1, Length[testdata]}]],
 {a, b, sigma}
]

As you can see, this returns a Function that takes a list as its first argument.

On a different note: when posting code, you can easily style it with Markdown to make it easier for someone else to copy-paste it.

[spasetto]

Thanks for your prompt reply! "Almost" everything working now! I think one issue is solved, but still, the code doesn't complete the example with predictiveDistribution and regressionPlot1D. Let me see if I manage with Markdown:

The object

obj = defineInferenceProblem[
  "Data" -> testdata,
  "LogLikelihoodFunction" -> expressionToFunction[
     LogLikelihood[NormalDistribution[0, sigma], 
     Table[
      testdata[[i, 2]] - (a testdata[[i, 1]] + b), {i, 1, 
       Length[testdata]}]],
     {a, b, sigma}
    ],
  "Parameters" -> {{a, -10, 10}, {b, -10, 10}, {sigma, 0.1, 10}},
  "PriorDistribution" -> {"LocationParameter", "LocationParameter",
    "ScaleParameter"}
  ]

and nestedSampling result

results = nestedSampling[obj]

are working pretty well. But,

predictions = 
  predictiveDistribution[takePosteriorFraction[results, 0.99], 
   Range[-2, 2, 0.1]];
regressionPlot1D[results, predictions]

is returning an empty box.

Vice versa, if I define the obj with as in your example

Clear[results,obj]

obj = defineInferenceProblem[
  "Data" -> testdata,
  "GeneratingDistribution" -> NormalDistribution[
    a x + b,
    sigma
    ],
  "Parameters" -> {{a, -10, 10}, {b, -10, 10}, {sigma, 0.1, 10}},
  "IndependentVariables" -> {x},
  "PriorDistribution" -> {"LocationParameter", "LocationParameter", 
    "ScaleParameter"}
  ]

...then

results = nestedSampling[obj]
predictions = 
  predictiveDistribution[takePosteriorFraction[results, 0.99], 
   Range[-2, 2, 0.1]];
regressionPlot1D[results, predictions]

is returning a pretty plot. Is there something wrong with predictiveDistribution or regressionPlot1D, or (most probably) did I miss a passage?

ThX in advance for your kind help!

[SjoerdSmitWolfram]

No problem. Unfortunately, the function predictiveDistribution needs to know the "GeneratingDistribution" property of your object. It's possible to do nested sampling without the "GeneratingDistribution" property, but finding the predictive distribution needs it. I just noticed that the warning message for this doesn't fire for this particular case. What you should see is:

Message[predictiveDistribution::MissGenDist]

Note that predictiveDistribution returns unevaluated in your example. You can still add the "GeneratingDistribution" to the inference object if you need predictiveDistribution to work. When you specify your own loglikelihood function, that's the one that will be used in the nested sampling algorithm and the "GeneratingDistribution" property is effectively ignored by defineInferenceProblem.

[spasetto]

Hi,

Thanks for your help with this. What do you mean for "You can still add the "GeneratingDistribution" to the inference object if you need predictiveDistribution to work"? I've tried

obj = defineInferenceProblem[
  "Data" -> testdata,
  "LogLikelihoodFunction" -> 
   expressionToFunction[
    LogLikelihood[NormalDistribution[0, sigma], 
     Table[
      testdata[[i, 2]] - (a testdata[[i, 1]] + b), {i, 1, 
       Length[testdata]}]], {a, b, sigma}],
  "Parameters" -> {{a, -10, 10}, {b, -10, 10}, {sigma, 0.1, 10}}, 
  "PriorDistribution" -> {"LocationParameter", "LocationParameter", 
    "ScaleParameter"}]

results = nestedSampling[obj]

but

 Append[results, 
 "GeneratingDistribution" -> NormalDistribution[a x + b, sigma]]

does not work to fix predictiveDistribution and regressionPlot1D.

predictions = predictiveDistribution[results, Range[-2, 2, 0.1]];
regressionPlot1D[results, predictions]

Still returns empty. Pity, because predictiveDistribution and regressionPlot1D seem so convenient!

ThX

[SjoerdSmitWolfram]

Append does not modify the variable results; that's not how WL works. Try:

results = Append[results, "GeneratingDistribution" -> NormalDistribution[a x + b, sigma]]

instead.

You can also just include "GeneratingDistribution" -> NormalDistribution[a x + b, sigma] in defineInferenceProblem.

[spasetto]

Hi, I'm really sorry. None of your options worked. Let me recap: for example you did

BlockRandom[SeedRandom[1];
 testdata = (# -> (# + 2 RandomReal[])) & /@ RandomReal[{-1, 1}, 20];
 ListPlot[List @@@ testdata]]

obj = defineInferenceProblem[
  "Data" -> testdata,
  (*"GeneratingDistribution"->NormalDistribution[a x+b,sigma],*)

  "LogLikelihoodFunction" -> 
   expressionToFunction[
    LogLikelihood[NormalDistribution[0, sigma], 
     Table[
      testdata[[i, 2]] - (a testdata[[i, 1]] + b), {i, 1, 
       Length[testdata]}]], {a, b, sigma}],
  "Parameters" -> {{a, -10, 10}, {b, -10, 10}, {sigma, 0.1, 10}}, 
  "PriorDistribution" -> {"LocationParameter", "LocationParameter", 
    "ScaleParameter"}]

results = nestedSampling[obj]
results = 
 Append[results, 
  "GeneratingDistribution" -> NormalDistribution[a x + b, sigma]]

predictions = predictiveDistribution[results, Range[-2, 2, 0.1]];
regressionPlot1D[results, predictions]

and you did get back the predictiveDistribution in regressionPlot1D plot (after 1 Listplot and 3 infer.objs)? (Mathematica 12.3.0.0, right?)

[ssmit1986]

Ah, I think I see why: it still needs the "IndependentVariables" -> {x} property to find the predictiveDistribution. It should work if you add that. I tried it on V12.3 as well.

[spasetto]

Hi, Sorry, I thought I closed already. It worked indeed! Thanks a lot! Best Stefano