If it helps, here is equivalent R code:

# Transmitter location
  a <- c(0.75, 0.8)

# Receiver locations
  n <- 4  # Number of recieviers
  r <- matrix(c(0.25, 0.25, 0.25, 0.75, 0.75, 0.25, 0.75, 0.75), nrow=4, byrow=TRUE)

# True squared distances to transmitter
  trueDist2 = (r[,1] - a[1])^2 + (r[,2] - a[2])^2
  # [1] 0.5525 0.2525 0.3025 0.0025

# Error parameter associated with measurement of squared distance
  truesigma <- 0.05

# Sample of observations of squared distance
  set.seed(12345)
  d2 <- rgamma(n, shape=trueDist2/truesigma^2, scale=truesigma^2) 
  d2 <- c(0.454823, 0.24621, 0.286231, 0.00226926)  # Squared distances to match Mathematica example
  
# Define log likelihood function
  logL <- function(parms, r, d2) {
    mx <- parms[1]
    my <- parms[2]
    sigma <- parms[3]
    temp <- ((r[, 1] - mx)^2 + (r[, 2] - my)^2)/sigma^2
    -sum(d2)/sigma^2 - sum(lgamma(temp)) - sum(log(sigma^2)*temp) + sum(log(d2)*(-1 + temp))
  }                
   
# Maximize log likelihood
  result <- optim(c(a[1], a[2], truesigma), logL, r=r, d2=d2, 
    control=list(fnscale=-1), hessian=TRUE)
  result

# Covariance matrix and standard errors for parameters
  cov <- -solve(result$hessian)
  se <- diag(cov)^0.5 
  # {0.01457950738463684, 0.01626917105787927, 0.01464057953398463} *)     

I think you’ll need to assign a distribution to the observed distance (or the square of the distance) that has a mean and variance dependent on the true distance. I don’t know the characteristics of the devices you’re using and offer the following approach only as a possible structure to follow.

I’m assuming from your description that from each receiver you only obtain a distance and NOT a direction.

Suppose the distribution of the square of the observed distance follows a Gamma distribution with the mean equal to the square of the true distance and the variance is proportional to the square of the true distance. (I’m assuming that smaller distances are more precisely estimated than large differences.)

If the true squared distance is $\mu^2$ and the variance is $\mu^2 \sigma^2$, then the associated Gamma distribution parameters are $\alpha=\mu^2/\sigma^2$ and $\beta=\sigma^2$.

Below is Mathematica code for generating data and finding the maximum likelihood estimates.

(* Transmitter location *)
a = {0.75, 0.8};

(* Receiver locations *)
r = {{0.25, 0.25}, {0.25, 0.75}, {0.75, 0.25}, {0.75, 0.75}};
n = Length[r] (* Number of receivers *)
trueDist2 = Table[Norm[r[[i]] - a]^2, {i, n}] (* True distances to transmitter *)
trueσ = 0.05;  (* Error parameter *)

(* Sample of observations *)
SeedRandom[12345];
d = Table[RandomVariate[GammaDistribution[trueDist2[[i]]/trueσ^2, trueσ^2]], {i, n}]
(* {0.5525,0.2525,0.3025,0.0025} *)

(* Log of the likelihood *)
logL = Sum[LogLikelihood[GammaDistribution[((r[[i, 1]] - μx)^2 + (r[[i, 2]] - μy)^2)/σ^2, σ^2], {d[[i]]}],
   {i, n}];

(* Maximum likelihood estimates *)
mle = FindMaximum[{logL, σ > 0}, {{μx, a[[1]]}, {μy, a[[2]]}, {σ, trueσ}}, MaxIterations -> 5000]
(* {11.571356964839643,{μx -> 0.7146660584863823, μy -> 0.7750478956280523, σ -> 0.04110551720224718}} *)

(* Estimate of covariance matrix *)
(cov = -Inverse[D[logL, {{μx, μy, σ}, 2}] /. mle[[2]]]) // MatrixForm
se = Diagonal[cov]^0.5  (* Standard errors of estimates of μx, μy, and σ *)
(* {0.01457950738463684, 0.01626917105787927, 0.01464057953398463} *)

I think you’ll need to assign a distribution to the observed distance (or the square of the distance) that has a mean and variance dependent on the true distance. I don’t know the characteristics of the devices you’re using and offer the following approach only as a possible structure to follow.

I’m assuming from your description that from each receiver you only obtain a distance and NOT a direction.

Suppose the distribution of the square of the observed distance follows a Gamma distribution with the mean equal to the square of the true distance and the variance is proportional to the square of the true distance. (I’m assuming that smaller distances are more precisely estimated than large differences.)

If the true squared distance is $\mu^2$ and the variance is $\mu^2 \sigma^2$, then the associated Gamma distribution parameters are $\alpha=\mu^2/\sigma^2$ and $\beta=\sigma^2$.

Below is Mathematica code for generating data and finding the maximum likelihood estimates.

(* Transmitter location *)
a = {0.75, 0.8};

(* Receiver locations *)
r = {{0.25, 0.25}, {0.25, 0.75}, {0.75, 0.25}, {0.75, 0.75}};
n = Length[r] (* Number of receivers *)
trueDist2 = Table[Norm[r[[i]] - a]^2, {i, n}] (* True squared distances to transmitter *)
trueσ = 0.05;  (* Error parameter *)

(* Sample of observations of squared distance *)
SeedRandom[12345];
d2 = Table[RandomVariate[GammaDistribution[trueDist2[[i]]/trueσ^2, trueσ^2]], {i, n}]
(* {0.5525,0.2525,0.3025,0.0025} *)

(* Log of the likelihood *)
logL = Sum[LogLikelihood[GammaDistribution[((r[[i, 1]] - μx)^2 + (r[[i, 2]] - μy)^2)/σ^2, σ^2], {d2[[i]]}],
   {i, n}];

(* Maximum likelihood estimates *)
mle = FindMaximum[{logL, σ > 0}, {{μx, a[[1]]}, {μy, a[[2]]}, {σ, trueσ}}, MaxIterations -> 5000]
(* {11.571356964839643,{μx -> 0.7146660584863823, μy -> 0.7750478956280523, σ -> 0.04110551720224718}} *)

(* Estimate of covariance matrix *)
(cov = -Inverse[D[logL, {{μx, μy, σ}, 2}] /. mle[[2]]]) // MatrixForm
se = Diagonal[cov]^0.5  (* Standard errors of estimates of μx, μy, and σ *)
(* {0.01457950738463684, 0.01626917105787927, 0.01464057953398463} *)