how to fit pdf of known form to data I have a set $X$ of 1000 data points.
I know the PDF has a certain form, but there are two constant parameters for which I need to derive values in order to bet fit the data.
Is there an established way to do this?
Off the top of my head (and, full disclosure, I'm a newbie, so this could be a terrible idea):


*

*Partition the data points into "bins".  (How many bins?  How to decide the bin ranges?)

*Use non-linear methods to derive constants that minimize the RMSE between the predicted # of data points per bin and the actual # of data points per bin.


Is that crazy?  Is there a better and/or simpler way?
Given I'm starting from scratch and know nothing about R and MATLAB, is one or the other particularly well-suited to solving this sort of problem?
I've written C code that used the gnu GSL library to do non-linear fitting, but surely there's a better way.
EDITED TO ADD:
If $P(x,m,s)$ is the PDF of a normal distribution with mean $m$ and stdev $s$ then my PDF is of the form
$P'(x,s,k) = (1 - 2k) P(x,0,s) + k  P(x,-1,s) + k  P(x,+1,s)$
where $0 \le k \le 0.5$. 
Basically there was an original series $X$ that was modified by shifting some of its elements by $+1$ and an equal number by $-1$.  All I know about the original series $X$ is that it had a normal distribution with mean $0$.
EDITED TO ADD:
Including an adapted version of user777's R snippet for posterity.  I kept running into the problem where the optimizer would choose very small values for sigma, which resulted in the likelihood of certain data points underflowing to zero, the log of which is -Inf, which caused the optimizer to barf.  To solve that I transform the parameters given to the function being optimized by mapping them onto a fixed range (using arctan). The optimizer can pass very large (or very small) values to the function and they'll be transformed into "reasonable" values.
n <- 1000
k <- 0.25

set.seed(0)

x <- rnorm(n, sd=0.25)

for(i in 1:floor(n*k*2)) 
{
    x[i] = x[i] + if (i %% 2) +1 else -1
}

plot(density(x))

klo = 0.0
khi = 0.5
slo = 0.0
shi = sd(x)  # sd of original series <= sd of modified series

at <- function(x, lo, hi) {
    return((atan(x) + pi/2 + lo)*(hi/pi))
}

ll <- function(params, data) {
    k <- at(params[1], klo, khi)
    s <- at(params[2], slo, shi)
    p <- (1 - 2 * k) * dnorm(data, mean=0, sd=s) + k * dnorm(data, mean=1, sd=s) + k * dnorm(data, mean=-1, sd=s)
    return(0 - sum(log(p)))
}

methods <- c("Nelder-Mead", "BFGS", "CG", "SANN")

for(m in methods) {
    result <- optim(c(0, 0), fn=ll, data=x, method=m)
    print(c(m, at(result$par[1], klo, khi), at(result$par[2], slo, shi)))
}

 A: The additional information you've added indicates that this is a mixture model (several PDFs composed together). These can also be estimated by MLE. Each data point has probability given by
$$
p(x|k,s)=(1-2k)\mathcal{N}(0,s)+k\mathcal{N}(1,s)+k\mathcal{N}(-1,s)
$$
But MLE maximizes the joint probability of all of these data points. So take the product, or, to mitigate the risk of underflow, take the sum of the logs, which is the same as taking the log of the product by the properties of logarithms.
So all you need to do is use a decent optimizer to find the maximum of the objective function. There's no guarantee that the maximum is unique in this context, so you actually have a global optimization problem. But finding each optima is cheap, so I'd recommend just using multistart with a fast optimizer.
N <- 1e5
k <- 1/10

set.seed(13)
x1 <- rnorm(N*(1-2*k), sd=3)
x2 <- rnorm(N*k, mean=1, sd=3)
x3 <- rnorm(N*k, mean=-1, sd=3)

x <- c(x1, x2, x3)

plot(density(x))

ll <- function(params, data){
    k <- params[1]
    sigma <- params[2]
    lik <- (1-2*k)*dnorm(data, mean=0, sd=sigma)+k*dnorm(data, mean=1, sd=sigma)+k*dnorm(data, mean=-1, sd=sigma)
    out <- sum(log(lik))
    return(out)
}

optim(runif(2, 0, 10), fn=ll, control=list(fnscale=-1), data=x, method="L-BFGS-B", lower=c(0,0), upper=c(0.5, Inf)) 

The results are ok. The estimate of the mixture coefficient $k$ is decidedly wrong, but the variance estimates are correct.
$par
[1] 0.1808578 2.9688096

$value
[1] -252721.3

$counts
function gradient 
      22       22 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

A: You can use Maximum Likelihood Estimation (MLE) to find the parameter of the distribution. MLE is typically done using Expectation Maximization (EM). In the following I explain what needs to be done to solve the problem, in more detail.
What needs to be done? Let's formalize the problem first.
Let $\mathcal{D} = \{x_1, \dots, x_n\}$ be the set of sample points; $n=1000$ in your case. Also, let $\theta$ denote the parameters of the target distribution; $\theta = (k, s)$ in your case. 
MLE objective is to find the parameters $\theta^*$ that maximize the likelihood function, which we denote by $L_\mathcal{D}$); that is:
$$
\theta^* = \arg\max_{\hspace{-0.75cm}\theta = (s, k)} L_\mathcal{D}(\theta) = \arg\max_{\hspace{-0.75cm}\theta = (s, k)} \log \left( L_\mathcal{D}(\theta) \right)
$$
The second equality holds because $\log$ is a monotonically increasing function. If we assume data samples are independent and identically distributed (i.i.d) we can write the likelihood function as:
$$
L_\mathcal{D}(\theta) = \prod_{i=1}^n P(x_i | \theta)
$$
Using this we can rewrite the MLE objective as follows:
$$
\theta^* = \arg\max_{\hspace{-0.75cm}\theta = (s, k)} \sum_{i=1}^n \log \left( P(x_i | \theta) \right)
\quad \quad (1)
$$
In your problem, the distribution $P(x | \theta)$ is a Gaussian Mixture Model with 3 mixture components. This can be defined as follows:
$$
P(x | \theta) = \sum_{z = 1}^3 P(x, z | \theta) = \sum_{z = 1}^3 P(z | \theta) P(x | \theta, z)\\ 
P(z | \theta) = 
\begin{cases}
k & \text{if $z=1$} \\
k & \text{if $z=2$} \\
2k-1 & \text{if $z=3$}
\end{cases}
, \quad P(x | \theta, z) = 
\begin{cases}
\mathcal{N}(0, s) & \text{if $z=1$} \\
\mathcal{N}(1, s) & \text{if $z=2$} \\
\mathcal{N}(-1, s) & \text{if $z=3$}
\end{cases}
$$
In summary, you need to solve the following optimization problem:
$$
\theta^* = \arg\max_{\hspace{-0.75cm}\theta = (s, k)} \sum_{i=1}^n 
\log \left( \sum_{z_i=1}^k P(z_i | \theta) P(x_i | \theta, z_i) \right)
\quad \quad (2)
$$
How can we solve the problem in practice?
(2) can be solved using Expectation Maximization (EM). EM is an iterative method which starts from an initial solution $\theta^{(0)}$ and updates the solution in each iteration. Let $\theta^{(t)}$ denote the solution of EM in the $t$-th iteration.
Since $\log$ is a concave function, we can use Jensen's Inequality to lower bound the objective function of (2) (which is the likelihood function) as follows:
$$
A_q(\theta) = \sum_{i=1}^n q_i(z_i) \log \left( P(x_i | \theta, z_i) \right)
\leq 
\sum_{i=1}^n 
\log \left( \sum_{z_i=1}^k P(z_i | \theta) P(x_i | \theta, z_i) \right)
$$
where $q = (q_1, \dots, q_n)$ and $\forall i, q_i$ is an arbitrary discrete distribution on the unobserved random variable $z_i$; that is, $q$ can be represented by an $n \times 3$ matrix with non-negative entries where the $i$-th row (denoted by $q_i$) sums to 1 (i.e. $\sum_{j=1}^3 q_{i,j} = 1$) and $q_i$ represents the distribution over 3 possible outcomes for $z_i \in \{1, 2, 3\}$.
Although this may not be obvious, but, it can be shown that for any given $\theta$, if we set $q_{i, j} = p(z_i = j | \theta, x_i)$ we have $A_q(\theta) = L_\mathcal{D}(\theta)$. The EM algorithm uses this property to ensure that the sequence of the solutions that it finds is increasing that is $\forall t > 0, L_\mathcal{D}(\theta^{(t)}) \geq L_\mathcal{D}(\theta^{(t-1)})$. The EM algorithm works as follows:


*

*initialize $\theta^{(0)}$; e.g. set it to $(0.25, 1)$ in your case

*set $t=0$

*fix $q_{i,j} = p(z_i = j | \theta^{(t)}, x_i)$

*update the solution $\theta^{(t+1)} = \arg\max_{\theta} A_q(\theta)$

*$t = t+1$

*if($\theta^{(t)} \neq \theta^{(t-1)}$) go to 3 else output $\theta^{(t-1)}$

