Survival estimation when death/censoring is probabilistic I am trying to estimate survival function, but in case where each event is censoring with probability $p_i$. (That is, I am never sure if the event is right-censoring or death, but I can estimate the probability of each.)
Typical Kaplan-Meier works only for binary right-censoring. Is there any straightforward generalization? 
 A: If we can assume that for each item/individual you have $ x_i $ (covariates), $ t_i $ (last observation) and $ p_i $ (probability of death rather than censoring), and that the $ p_i $ are conditionally independent given $ x_i $ and $ t_i $: I don't see a problem with a simple modification to the Kaplan-Meier product limit formula. 
If the standard formula is:
$$
\hat S (t) = \prod_{t_i<t}{\frac{n_i-d_i}{n_i}}
$$
Where $ n_i $ is the number of individuals at risk at time $ t_i $ and $ d_i $ is the number dying at time $ t_i $.
Then you should be able to use
$$
\hat S (t) = \prod_{t_i <t}{\frac{n_i-\sum_{t_j=t_i}{p_j}}{n_i}}
$$
Where $ n_i $ is now the expected number of individuals at risk. (Note that I have actually used $ p_i $ as probability of death rather than censoring for notational convenience.)
For other analyses (e.g., Cox proportional hazards) you may find it more difficult to construct the likelihood function in a computationally convenient or tractable form. What you should be able to do if nothing else is sample $ D_i \sim Bernoulli (p_i) $ and perform your analyses multiple times in a Monte Carlo manner. 

EDIT
I have now simulated this in R and I am happy with it. Briefly I created the following data generating process:
\begin{align}
T_{death} &\sim Weibull(\alpha_{death},\beta_{death}) \\
T_{cens}  &\sim Weibull(\alpha_{cens}, \beta_{cens} ) \\
T &= \min\{T_{death},T_{cens}\} \\
\lambda_{death}(t) &= \frac{\alpha_{death}}{\beta_{death}}\left(\frac{t}{\beta_{death}}\right)^{\alpha_{death}-1} \\
\lambda_{cens}(t) &= \frac{\alpha_{cens}}{\beta_{cens}}\left(\frac{t}{\beta_{cens}}\right)^{\alpha_{cens}-1} \\
p(t) &= \frac{\lambda_{death}(t)}{\lambda_{death}(t)+\lambda_{cens}(t)}
\end{align}
Using the estimator as I suggested above gives the following good match to the survival curve for $T_{death}$ (i.e., the censoring process has not interfered with estimation).

The code for it is:
library(ggplot2)

# Data generating function parameters
a_death = 1.5
b_death = 10.0
a_cens = 1.0
b_cens = 7.5
n = 1000

# Generate data
t <- c(rweibull(n*mix_deathcens, a_death, b_death), rweibull(n*(1-mix_deathcens), a_cens, b_cens))
t_death <- rweibull(n, a_death, b_death)
t_cens  <- rweibull(n, a_cens,  b_cens )
t <- pmin(t_death, t_cens)
g <- (t_death < t_cens) * 1

# Assign probabilities based on hazard rates
hr_death <- (a_death/b_death) * (t/b_death)^(a_death-1)
hr_cens  <- (a_cens /b_cens ) * (t/b_cens )^(a_cens -1)
p <- hr_death/(hr_death+hr_cens)

# PRODUCT LIMIT
ord <- order(t)
df <- data.frame(t=t[ord], p=p[ord], S=rep(NA,n))
S_i <- 1
n_i <- n
for (i in 1:length(df$t)) {
  n_i <- n_i - (1 - df$p[i])
  df$S[i] <- S_i * (n_i - df$p[i])/n_i
  S_i <- df$S[i]
  n_i <- n_i - df$p[i]
}

# Calculate desired survival curve
df$S_death <- pweibull(df$t, a_death, b_death, lower.tail = FALSE)

# Plot
ggplot(df) + geom_step(aes(x = t, y = S)) + geom_line(aes(x = t, y = S_death)) + xlab("Time") + ylab("Survival")

